{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Recursive least squares\n",
    "\n",
    "Recursive least squares is an expanding window version of ordinary least squares. In addition to availability of regression coefficients computed recursively, the recursively computed residuals the construction of statistics to investigate parameter instability.\n",
    "\n",
    "The `RecursiveLS` class allows computation of recursive residuals and computes CUSUM and CUSUM of squares statistics. Plotting these statistics along with reference lines denoting statistically significant deviations from the null hypothesis of stable parameters allows an easy visual indication of parameter stability.\n",
    "\n",
    "Finally, the `RecursiveLS` model allows imposing linear restrictions on the parameter vectors, and can be constructed using the formula interface."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "from pandas_datareader.data import DataReader\n",
    "\n",
    "np.set_printoptions(suppress=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 1: Copper\n",
    "\n",
    "We first consider parameter stability in the copper dataset (description below)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(sm.datasets.copper.DESCRLONG)\n",
    "\n",
    "dta = sm.datasets.copper.load_pandas().data\n",
    "dta.index = pd.date_range(\"1951-01-01\", \"1975-01-01\", freq=\"YS\")\n",
    "endog = dta[\"WORLDCONSUMPTION\"]\n",
    "\n",
    "# To the regressors in the dataset, we add a column of ones for an intercept\n",
    "exog = sm.add_constant(\n",
    "    dta[[\"COPPERPRICE\", \"INCOMEINDEX\", \"ALUMPRICE\", \"INVENTORYINDEX\"]]\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, construct and fit the model, and print a summary. Although the `RLS` model computes the regression parameters recursively, so there are as many estimates as there are datapoints, the summary table only presents the regression parameters estimated on the entire sample; except for small effects from initialization of the recursions, these estimates are equivalent to OLS estimates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod = sm.RecursiveLS(endog, exog)\n",
    "res = mod.fit()\n",
    "\n",
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The recursive coefficients are available in the `recursive_coefficients` attribute. Alternatively, plots can generated using the `plot_recursive_coefficient` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(res.recursive_coefficients.filtered[0])\n",
    "res.plot_recursive_coefficient(range(mod.k_exog), alpha=None, figsize=(10, 6))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The CUSUM statistic is available in the `cusum` attribute, but usually it is more convenient to visually check for parameter stability using the `plot_cusum` method. In the plot below, the CUSUM statistic does not move outside of the 5% significance bands, so we fail to reject the null hypothesis of stable parameters at the 5% level."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(res.cusum)\n",
    "fig = res.plot_cusum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another related statistic is the CUSUM of squares. It is available in the `cusum_squares` attribute, but it is similarly more convenient to check it visually, using the `plot_cusum_squares` method. In the plot below, the CUSUM of squares statistic does not move outside of the 5% significance bands, so we fail to reject the null hypothesis of stable parameters at the 5% level."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res.plot_cusum_squares()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 2: Quantity theory of money\n",
    "\n",
    "The quantity theory of money suggests that \"a given change in the rate of change in the quantity of money induces ... an equal change in the rate of price inflation\" (Lucas, 1980). Following Lucas, we examine the relationship between double-sided exponentially weighted moving averages of money growth and CPI inflation. Although Lucas found the relationship between these variables to be stable, more recently it appears that the relationship is unstable; see e.g. Sargent and Surico (2010)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start = \"1959-12-01\"\n",
    "end = \"2015-01-01\"\n",
    "m2 = DataReader(\"M2SL\", \"fred\", start=start, end=end)\n",
    "cpi = DataReader(\"CPIAUCSL\", \"fred\", start=start, end=end)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ewma(series, beta, n_window):\n",
    "    nobs = len(series)\n",
    "    scalar = (1 - beta) / (1 + beta)\n",
    "    ma = []\n",
    "    k = np.arange(n_window, 0, -1)\n",
    "    weights = np.r_[beta**k, 1, beta ** k[::-1]]\n",
    "    for t in range(n_window, nobs - n_window):\n",
    "        window = series.iloc[t - n_window : t + n_window + 1].values\n",
    "        ma.append(scalar * np.sum(weights * window))\n",
    "    return pd.Series(ma, name=series.name, index=series.iloc[n_window:-n_window].index)\n",
    "\n",
    "\n",
    "m2_ewma = ewma(np.log(m2[\"M2SL\"].resample(\"QS\").mean()).diff().iloc[1:], 0.95, 10 * 4)\n",
    "cpi_ewma = ewma(\n",
    "    np.log(cpi[\"CPIAUCSL\"].resample(\"QS\").mean()).diff().iloc[1:], 0.95, 10 * 4\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After constructing the moving averages using the $\\beta = 0.95$ filter of Lucas (with a window of 10 years on either side), we plot each of the series below. Although they appear to move together prior for part of the sample, after 1990 they appear to diverge."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(13, 3))\n",
    "\n",
    "ax.plot(m2_ewma, label=\"M2 Growth (EWMA)\")\n",
    "ax.plot(cpi_ewma, label=\"CPI Inflation (EWMA)\")\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "endog = cpi_ewma\n",
    "exog = sm.add_constant(m2_ewma)\n",
    "exog.columns = [\"const\", \"M2\"]\n",
    "\n",
    "mod = sm.RecursiveLS(endog, exog)\n",
    "res = mod.fit()\n",
    "\n",
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res.plot_recursive_coefficient(1, alpha=None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The CUSUM plot now shows substantial deviation at the 5% level, suggesting a rejection of the null hypothesis of parameter stability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res.plot_cusum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, the CUSUM of squares shows substantial deviation at the 5% level, also suggesting a rejection of the null hypothesis of parameter stability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res.plot_cusum_squares()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 3: Linear restrictions and formulas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear restrictions\n",
    "\n",
    "It is not hard to implement linear restrictions, using the `constraints` parameter in constructing the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "endog = dta[\"WORLDCONSUMPTION\"]\n",
    "exog = sm.add_constant(\n",
    "    dta[[\"COPPERPRICE\", \"INCOMEINDEX\", \"ALUMPRICE\", \"INVENTORYINDEX\"]]\n",
    ")\n",
    "\n",
    "mod = sm.RecursiveLS(endog, exog, constraints=\"COPPERPRICE = ALUMPRICE\")\n",
    "res = mod.fit()\n",
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Formula\n",
    "\n",
    "One could fit the same model using the class method `from_formula`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod = sm.RecursiveLS.from_formula(\n",
    "    \"WORLDCONSUMPTION ~ COPPERPRICE + INCOMEINDEX + ALUMPRICE + INVENTORYINDEX\",\n",
    "    dta,\n",
    "    constraints=\"COPPERPRICE = ALUMPRICE\",\n",
    ")\n",
    "res = mod.fit()\n",
    "print(res.summary())"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
